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ABSTRACT 

In this first paper of this series, we present a new approach for studying the chemodynamical evolu- 
tion in disc galaxies, which consists of fusing disc chemical evolution models with compatible numerical 
simulations of galactic discs. Such a method avoids known star formation and chemical enrichment 
problems encountered in simulations. Here we focus on the Milky Way, by using a detailed thin-disc 
chemical evolution model (matching local observables, which are weakly affected by radial migration) 
and a simulation in the cosmological context, with dynamical properties close to those of our Galaxy. 
We examine in detail the interplay between in-situ chemical enrichment and radial migration, and 
their impact on key observables in the solar neighborhood, e.g., the age-metallicity- velocity relation, 
the metallicity distribution, and gradients in the radial and vertical directions. We show that, due to 
radial migration from mergers at high redshift and the central bar at later times, a sizable fraction of 
old metal-poor, high-[a/Fe] stars can reach the solar vicinity. This naturally accounts for a number of 
observations related to both the thin and thick discs, despite the fact that we use thin-disc chemistry 
only. Although significant radial mixing is present, a slope in the age-metallicity relation is preserved, 
with a scatter compatible with recent observational work. While we find a smooth density distribution 
in the [0/Fe]-[Fe/H]-plane, we can recover the observed discontinuity by selecting particles according 
to kinematical criteria used in high-resolution samples to define the thin and thick discs. We show 
that in the absence of early-on massive mergers the vertical velocity dispersion of the oldest stars is 
underestimated by a factor of ~ 2 compared to observations. We predict that the most likely birth 
place for the Sun is in the range 4.4 < r < 7.7 kpc (for a current location at r = 8 kpc), with the 
highest probability *~5.6 kpc, followed by ~7 kpc. Finally, we offer a new, unifying model for the 
Milky Way thick disc, where both mergers and radial migration play a role at different stages of the 
disc evolution. 

Subject headings: ISM: abundances — Galaxy: abundances — galaxies: ISM — Galaxy: evolution — 
galaxies: kinematics and dynamics — galaxies: structure 



1. INTRODUCTION 

Crucial information regarding the dominant mecha- 
nism responsible for the formation of the Milky Way 
(MW) disc and other Galactic components, is encoded 
in the chemical and dynamical properties of its stars. 
This conviction has led to an unprecedented observa- 
tional efforts in the past decades, aimed at mapping 
the chemistry and kinematics of a large number of stars 
in the MW. Until the end of 2003, most of the infor- 
mation was confined to small local samples, for which 
high-resolu tion s pectro scopic data was o btained (e.g., 
Fuhrmannl 119981 120081 within 25 pc and IBensbv et all 
20031 within 1 00 pc). In 2004 the G eneva Copenhagen 
Survey (CGS iNordstrom et all 120041 ) obtained the first 
large spectro-photometric sample of around 16 000 stars, 
as part of a Hipparcos follow-up campaign (hence, also 
confined to 100 pc from the Sun). More recently, opti- 
cal spectroscopic low-resol ution surveys, such as SEGUE 
(Yan nv and Rockosil[200l and RAVE (jSteinmetz et al.1 
2006), have extended the studied volume to distances of 
a few kpc from the Sun (with the majority of stars in 
the distance range 0.5-3 kpc), and increased the num- 
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bers of stars with chemo-kinematical information by 
more than an order of magnitude (>200 000 spectra for 
SEGU E and >500 000 spectra for RAVE, see ISteinmetzl 
120121 ). This effort will be soon complemented by high- 
resolution spectroscopic survey s both in the optical - 
HERMES (Fr eeman et alll20"Tol) and in the near- infrared 
-APOGEE (jAllende Prieto et alj|200l iMaiewski et~aTl 
|2010T) . APOGEE aims at measuring chemo-kinematical 
properties of around 10 5 stars close to the galactic plane, 
thus complementing SEGUE and RAVE (which exclude 
most stars below ^200 pc), whereas HERMES aims at 
obtaining chemical informa tion for around 10 6 stars. In 
the near future, 4MOST (jde Jong et all 120121 ), proba- 
bly the most ambitious project, aims at sampling even 
larger volumes by obtaining chemo-kinematical proper- 
ties of many millions of stars (taking full advantage of the 
Gaia astrometric results) . The common aim of the huge 
observational campaigns briefly summarized above is to 
constrain the MW assembly history - one of the main 
goals of the newly emerged field, Galactic Archaeology. 

The underlying principle of Galactic Archaeology is 
that the chemical elements synthesized inside stars, and 
later ejected back into the interstellar medium (ISM), 
are incorporated into new generations of stars. As differ- 
ent elements are released to the ISM by stars of differ- 
ent masses and, therefore, on different timescales, stellar 
abundance ratios provide a cosmic clock, capable of elic- 
iting the past history of star formation and gas accretion 
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of a galaxjQ. One of the most widely used "chemical- 
clocks" is the [a/Fe] ratic0. 

1.1. Galactic Archeology and radial migration 

The power of Galactic Archaeology has been threat- 
ened both by observational and theoretical results, show- 
ing that stars most probably move away from their 
birthplaces, i.e, migrate radially. Observational signa- 
tures of this radial migration (or mixing) have been re- 
ported in the liter ature since th e 1970 's, with the pi- 
oneering works by iGrenonl (|1972t 119891 ) . Grcnon iden- 
tified an old population of super-metal-rich stars (here- 
after SMR) , presently at the Solar vicinity, but with kine- 
matics and abundance properties ind icative of an origin 
in t he inner Galactic di sc (see also iCastro et al.l 119971 
an d iTrevisan et aL] 1201 lh . These results were extended 
bv iHavwoodl (|2008[ ). who showed, by re-analyzing the 
GCS data, that the low- and high-metallicity tails of 
the thin disc are populated by objects, whose orbital 
properties suggest origin in the outer and inner Galac- 
tic disc, respectively. In particular, the so-called SMR 
stars show metallicities which exceed the present day 
ISM and th ose of young stars at t he solar vicinity. As 
discussed b y Chiap pini et al.l (|2003l ) (see also Table 5 by 
lAsplund et al.ll2009l ). the metallicitv at the solar vicinity 
is not expected to increase much since the Sun's forma- 
tion, i.e., in the last ~4 Gyr, due to the rather inefficient 
star formation rate (SFR) at the solar radius during this 
period, combined with co ntinuous gas infal l into the disc. 
Hence, as summarized in iChiappinil (|2009[ ). pure chemi- 
cal evolution models for the MW thin disc cannot explain 
stars more metal rich than ^0.2 dex and radial migration 
has to be invoked. 

N-body simulations hav e also long shown th at radial 
migration is unavoidable. iRaboud et al.l (|1998P l studied 
numerical sim ulations aimed at explaining the results by 
IGrenonl (fl989l) of a mean positive U-motion (where U is 
the Galactocentric radial velocity component of stars), 
which the authors interpreted as metal-rich stars from 
the inner galaxy, wandering i n the solar neighborhood. 
However. IRaboud et al.l (|1998| ) interpreted their findings 
as stars on hot bar orbits, not recognizing that perma- 
nent changes to the stellar a ngular momenta are possible . 
It was not until the work by Sellwood and Binnevl (2002) 
that radial migration was established as an important 
process affecting the entire disc, where stars shift guiding 
radii due to interaction with transient spi ral structure. 
Modern, h igh-resolution simulations (e.g., iRoskar et al.l 
120081 1201 ID have confirmed this finding, but left the role 
of the Galactic bar unexplored. 

A different radial migra ti on me chanis m was proposed 
bv IMinchev and Famaevl (|2010( ) and Mi nchev et al.1 
(|2011af l. resulting from the non-linear coupling between 
the bar and spiral waves, or spir a ls of different multiplic - 
ity (|Minchev and Quillenl 120061: IMinchev et~aH l2012al ). 
These works, al ong with studies of diffusion coefficient s 
in barred discs (|Brunetti et al.ll2Wl [Shevchenkdl20rl . 
predict a variation in migration efficiency with time and 

3 In most cases, the stellar surface abundances reflect the com- 
position of the interstellar medium at the time of their birth; this 
is the reason why stars can be seen as fossil records of the Galaxy 
evolution. 

4 Here we use the notation in brackets to indicate abundances 
relative to the Sun, i.e., [X/Y] = log(X/Y)- log(X/Y) . 



disc radius, establishing that the dynamical influence of 
the bar plays an integral part of MW d isc modeling. Re- 
cently, iComparetta and Quillenl (|2012l ) showed that ra- 
dial migration can result from short-lived density peaks 
arising from interference of spiral density waves, even 
if patterns are long-lived. Aside from internal structure, 
perturbations due to minor mergers have al so been shown 
to be effective at mix ing the outer discs (|Quillen et al.l 
20091 iBird etaL|[20T2T) . but also can, indirectly, affect the 
entire disc by inducing (reinforcing) spiral and bar insta- 
bilities. Considering the established presence of a central 
bar, spiral structure and evidence for merger activity in 
the MW, it is clear that all of the above mentioned ra- 
dial migration mechanisms would have an effect on the 
Galactic disc. 

In summary, a number of both observational and theo- 
retical results challenge the power of Galactic Archaeol- 
ogy. Therefore, the only possible way to advance in this 
field is the development of chemodynamical models tai- 
lored to the MW, in the cosmological framework. Only 
then, a meaningful comparison with the large amounts of 
current and forthcoming observational data (as summa- 
rized in the beginning of this Sec), can be carried out. 
This is the main goal of the present work, namely, to de- 
velop a chemodynamical model for the MW, to be able 
to quantify the importance of radial mixing throughout 
the evolution of our Galaxy. 

1.2. Difficulties with fully self- consistent simulations 

Producing disc-dominated galaxies has tradition- 
ally been challenging for cosmological models. In 
early simulations, extreme angular momentum loss 
during mergers gave birth to galaxies with overly- 
conce ntrated mass distribut i ons and massive bulges 
(e.g.. iNavarro and Benzl[l991l : iNavarro and White! Il994t 
lAbadi et all 1200.1 . Although an increase in res- 
olution and a better modeling of star formation 
and feedback have allowed recent simulations to pro- 
duce MW-mass gal a xies with reduced bulge fractions 
(lAgertz et al.l 120111 : iGuedes et all 120111 : iMartig et ail 
12012ft . none of these simulations include chemical evo- 
lution. Galaxy formation simulations including some 
treatment of ch emical evolution have been performed by 
several groups dRaiteri et ail 1199 6: Mos coni et al.l 120011: 
Lia et all l2002t Kawata and Gibson! I2003t~ Kobavashi| 



2004; Scannapicco et al. 2005; Martinez-Serrano et al. 
20081: iQppenheimer and Davel2008HWiersma et al.ll200£ 
Few et al.l I2012T ) . However, although the results are en- 



couraging, and global observed trends seem to be re- 
produced, such as th e mass-metallicity relation (e.g., 
Kobayashi et al. 20Q3) or the metallicity tr ends between 
the d ifferent galactic components (e.g., iTissera et al.1 
I2012T ) . it is still a challenge for these simulations to 
reproduce the properties of the MW (e.g ., the typical 
meta llicities of the different components - ITissera et al.l 
l2012f) . Additionally, the fr action of low metallicity 
stars are often overestima ted (|Kobavashi and Nakasatol 
120111 : iCalura et al.l l2012h . and reproducing the posi- 
tion of thin- and thick-disc stars in the [0/Fe l-[Fe/Hl 
plane has proved challenging (|Brook et al.ll2012T ) . While 
such issues could be due to unresolved metal mixing 
(|Wiersma et al.|[2009D , it is also worth noticing that none 
of the above-mentioned simulations reproduces simulta- 
neously the mass, the morphology and star formation 
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history (SFH) of the MW. 

This situation has led us to look for a novel way to 
approach this complex problem. We will show that this 
approach works encouragingly well, explaining not only 
current observations, but also leading to a more clear 
picture regarding the nature of the MW thick disc. 

1.3. Thick discs formation scenarios 

The large uncertainties in important observational con- 
straints in the MW, such as the age-velocity-metallicity 
relation, abundance gradients and their evolution, to- 
gether with the inherent complexity of the topic of 
Galaxy Assembly, have led to different scenarios to be 
proposed for the formation of the thick disc. 

One possibility is that thick discs were born thick at 
high redshift from the internal gravit ational instabili- 
ties in gas-rich, turbulen t, clumpy discs ([Bournaud et al.l 
2009; iForbes et al.l20H or in the turb ulent phase associ- 
ated with numerous gas-rich mergers ()Brook et al.ll2004 
120051) . They could also h ave been created through accre- 
tion of galaxy satellites (jMeza et al.l 120051 : lAbadi et al.l 
2003), where thick disc stars then have an extragalactic 
origin. 

Another possibility is that thick discs are created 
through the heating of preexisting thin discs. This can 
either happen fast, on a Gyr timescale, in an early vi- 
olent epoch, or as a continuous process throughout the 
galaxy lifetime. In the first case, the thick disc would 
appear as a clear distinct component in chemistry and 
phase-space, while in the latter case it would rather be 
seen as a gradual, continuous transition. 

Preexisting thin discs can be he ated fast through 
the effect of multiple minor mergers (|Quinn et al.l 119931 : 
iVillalobos and Helmil l2008t iDi Matteo et all 120111 ). the 
rate of which decreases with decreasing redshift. Evi- 
dence for merger encounters can b e seen in structure in 
phase-space of MW disc stars fe.g.. IMinchev et aT1l2009l : 
iGomez et alj|2012dlal). which could last for as long as 
~4 Gvr ([Gomez et al.ll2012bf ). 

Finally, a recently proposed mechanism for the forma- 
tion of thick discs is radial migration, which we discuss 
next. 



1.4. Radial migration and thick discs 

In the past several years there has been a growing con- 
viction that radial migration, driven by transient spirals, 
can be responsible for the formation of thick discs, by 
bringing out high-velocity-dispersion stars from the in- 
ner disc and the bulge (no need for mergers). This sce- 
nario was used, for ex a mple, in the analytical model of 
Schonrich and Binney (2009), where the authors man- 
aged to explain the MW thin- and thick-disc characteris- 
tics without the need of mergers or any discrete heating 
processes. Similarly, the increa se of disc thickness with 
time found in the simulation by iRoskar et al. _( 2008) has 
been a ttrib uted to migration in th e works bv lSales et all 
(2009]) and lLoebman et all (poTTh . 

However, how exactly radial migration affects disc 
thickening in dynamical models had not been demon- 
strated (but o n ly assumed), u ntil the recent work by 
iMinchev etUI (|2011bl l2012bl ). The latter authors 
showed unambiguously that radial migration driven by 



seculai0 (or internal) evolution has minor effect on disc 
thickening, mostly beyond 3 disc scale-lengths, thus, re- 
sulting in a flared disc0. This is due to the conservation 
of vertical action of migrators, as opposed to the incor- 
rect assumption of vertical energy conservation. While 
not building up a kinematical thick disc (since the effect 
of outward and inwar d migrators mostly cancels out), 
IMinchev et al.l (|2012bD noted that radial migration does 
contribute to a chemical thick disc, in the sense that 
outward migrators are preferentially deposited at higher 
distances above the galactic plane, with the converse ef- 
fect for inward migrators. It can be expected that mi- 
gration can thicken the disc, provided that stars were 
"preheated", e.g., either by mergers or being born hotQ 
(both of these are expected at high redshift). 

Once again, to advance our understanding of the MW 
disc formation and evolution, we need to make quantita- 
tive estimates of the importance of radial mixing, guided 
by observational constraints. This implies the need for 
chemodynamical models making predictions specifically 
for the solar vicinity, where most of the current obser- 
vational constraints are found. In Sec. [2] we describe 
the disc simulation we adopt in the present work and in 
Sec. [3] we describe our input chemistry. Sec. [4] is devoted 
to our new approach. Our results are shown in Sec. [SJ 
while a new explanation for the origin of the thick disc 
can be found in Sec. [Jj Here we concentrate on results 
for the solar vicinity, while results for the whole disc are 
shown in Paper II of this series. Conclusions are drawn 
in Sec. E 

2. A LATE-TYPE DISC GALAXY SIMULATION IN THE 
COSMOLOGICAL CONTEXT 

To properly model the MW, it is crucial to be con- 
sistent with some observational constraints at redshift 
z = , for example, a flat rotation curve (e.g., iMerrifieldl 
119921 ). a small bulge (e.g., iBinnev et al.lll997fl. a central 
bar o f an intermediate size (e.g. JBabusiaux and Gilmorel 
1200511 gas to tota l disc mass ratio of ~ 0.14 at the so- 
lar vicinity (e.g., iChiappini et al.l [20011 and references 
therein), and local disc velocity dis persions close to the 
observed ones (e.g. JLee et afll201ll) . 

It is clear that cosmological simulations would be the 
natural framework for a state-of-the-art chemodynami- 
cal study of the Milky Way. Unfortunately, as described 
in Sec. 11.21 a number of star formation and chemical en- 
richment problems still exist in fully self-consistent sim- 
ulations. We, therefore, resort to (in our view) the next 
best thing - a high-resolution simulation in the cosmo- 
logical context coupled with a pure chemical evolution 
model (Sec. [3]), as described in detail in Sec. HJbelow. 

The simulation used in this work is p art of a suite 
of numerical experiments first presented bv lMartig et all 

5 The expression "secular evolution" in this case is a misnomer, 
since in addition to "internal", it also has the meaning of "slow". 
In fact, radial migration happens on a radial epicyclic frequency 
time-scale, making it a very fast process. We, thus, simply use the 
term "internal" for the rest of the paper. 

6 IMin chev et al. (2012b) showed that migrators shrink (verti- 
cally cool) the disc inside the bar's CR, further contributing to 
disc flaring. This is related to the predominance of inward migra- 
tors at such small radii with a mean action smaller than that of 
the local population. 

7 Note that mig ration efficiency drops with velocity dispersion 
HSolwav et al.ll20T2 ). 
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Fig. 1. — First row: The left panel shows the rotational velocity (dashed curve) and circular velocity (solid curve) at the final simulation 
time. The middle panel presents the m = 2 Fourier amplitudes, A2/A0, as a function of radius estimated from the stellar density. Curves 
of different colors present the time evolution of A2/A0. To see better the evolution of the bar strength with time, in the right panel we 
show the amplitude averaged over the bar maximum. Second row: Face-on density maps of the stellar component for different times, as 
indicated. Third row: The corresponding edge-on view. Contour spacing is logarithmic Fourth row: Changes in angular momentum, 
AZ/, as a function of radius, estimated in a time window of 0.52 Gyr, centered on the times of the snapshots shown above. Both axes 
are divided by the circular velocity, thus units are kpc (galactic radius). Strong variations are seen with cosmic time due to satellite 
perturbations and increase in bar strength. 

scale the disc radius by a factor of 1.67. This places the 
bar's corotation resonance (CR) and 2:1 outer Lindblad 
resonance (OLR) at ~ 4.7 and ~ 7.5 kpc, respectively, 
consistent with a number o f stu died (e.g . . lDehnenll2000l 
iMinchev et ail [20071 [20101 see iGerhardl 120111 for a re- 
view). At the same time the disc scale- length, measured 
from particles of all ages in the range 3 < r < 15 kpc, 
beco mes ~ 3 kpc, in close agreement with observations 
(e.g. JGerhard1[200l . After this change^, our simulated 
disc satisfies the criteria outlined at the beginning of this 
Section, required for any dynamical study of the MW, in 
the following: 

(i) it has an approximately flat rotation curve with 
a circular velocity V c ~ 210 km/s at 8 kpc (slightly 
lower than the MW), shown in the first row, left panel 
of Fig. [TJ where we have corrected for asymmetric drift 



( 2012) , where the authors studied the evolution of 33 sim- 
ulated galaxies from z = 5 to z = using the zoom-in 
technique described in IMartig et al.l (|2009l ) . This tech- 
nique consists of extracting merger and accretion histo- 
ries (and geometry) for a given halo in a A-CDM cosmo- 
logical simulation, and then re-simulating these histories 
at much higher resolution (150 pc spatial, and 10 4-5 M 
mass resolut i on). The interested reader is referred to 
IMartig et all ([20121 ) for more information on the simula- 
tion method. 

The galaxy we have chosen has a number of properties 
consistent with the MW, including a central bar, as we 
describe below. We would like to explore the disc evo- 
lution for a time period of about 11 Gyr, which is close 
to the age of the oldest disc stars in the MW. However, 
none of the above galaxies start forming discs earlier than 
~ 9— 10 Gyr before z — 0. To remedy this, we extend our 
simulation for two additional Gyr, allowing for 11.2 Gyr 
of self-consistent evolution. 

To match the MW in terms of dynamics, we down- 



8 Variations in our results with different rescaling are discussed 
in the Appendix. 

9 Note that we rescale the disc radius so as to reflect the position 
of the bar's resonances at the final time. 
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as described in, e.g.. iBinnev and Tremainel (|2008T ) . 

(ii) the bulge is relatively small, with a bu lge-to-total 
ratio of ~l/5 (as measured with GALFIT - iPeng et al.l 
2002 - on a mock i-band im age, see bottom panel of 
Fig. 27 bv lMartig et~aI|[20Tl . 

(iii) it contains an intermediate size bar at the final 
simulation time, which develops early on and grows in 
strength during the disc evolution (see middle and right 
panels of Fig. Q] and discussion below). 

(iv) the disc grows self-consistently as the result of cos- 
mological gas accretion from filaments and (a small num- 
ber of) early-on gas-rich mergers, as well as merger de- 
bris, with a last significant merger concluding ^9 — 8 Gyr 
ago. 

(v) the disc gas-to-total mass ratio at the final time 
is ^0.12, consistent with the estimate of ~ 0.14 for the 
solar vicinity (Fig. [5J top right panel). 

(vi) the radial and vertical velocity dispersions at r w 
8 kpc are ~ 40 and ~ 20 km/s, in very good agreement 
with observations (see Sec. 15.51 and Fig. [6j) . 

We define t = as the time after the central bulge 
(spheroidal component) has formed. Only about 2% of 
these stars contribute to the stellar density at the final 
simulation time in the region 7 < r < 9 kpc, \z\ < 3 kpc, 
which we investigate in this paper. Therefore, for this 
work we ignore those old stars and focus on the particles 
forming mostly from gas infall at later times. 

In the first row, middle panel of Fig. [1] we plot the 
Fourier amplitude, A m j Aq, of the density of stars as a 
function of radius, where Aq is the axisymmetric com- 
ponent and to is the multiplicity of the pattern; here we 
only show the to = 2 component, A2/A0. Different col- 
ors indicate the evolution of the A2 radial profile in the 
time period specified by the color bar. The bar is seen 
to extend to ~ 3 — 4 kpc, where deviations from zero be- 
yond that radius are due to spiral structure. The brown 
curve reaching A2/A0 ~ 0.9 at r « 4.5 kpc results from a 
merger-induced two-armed spiral. To see better the evo- 
lution of the bar strength with time, in the right panel 
we show the amplitude averaged over 1 kpc at the bar 
maximum. 

The second and third rows of Fig. Q] show face-on and 
edge-on stellar density contours at different times of the 
disc evolution, as indicated in each panel. The redistri- 
bution of stellar angular momentum, L, in the disc as a 
function of time is shown in the fourth row, where AL is 
the change in the specific angular momentum as a func- 
tion of radius estimated in a time period At = 520 Myr, 
centered on the time of each snapshot above. Both axes 
are divided by the rotational velocity at each radius, 
therefore L is approximately equal to the initial radius 
(at the beginning of each time interval) and AL gives the 
distance by which guiding radii change during the time 
interval At. The dotted-red and solid-blue vertical lines 
indicate the positions of the bar's CR and OLR. Note 
that due to the bar's slowing down, these resonances are 
shifted outwards in the disc with time. 

After the initial bulge formation, the largest merger 
has a 1:5 stellar-mass ratio and an initially prograde or- 
bit, plunging through the center later and dissolving in 
~ 1 Gyr (in the time period 1.5 < t < 2.5 Gyr, first 
column of Fig. [T]). Due to its in-plane orbit (inclina- 
tion < 45°), this merger event results in accelerated disc 
growth by triggering strong spiral structure in addition 



to its tidal perturbation ()Quillen et al.l [2009). One can 
see the drastic effect on the changes in angular momen- 
tum, AL, in the fourth row, right panel of Fig. [1] More 
discussion about the importance of this merger can be 
found in Sections 16.31 and 19.21 A number of less violent 
events are present at that early epoch, with their fre- 
quency decreasing with time. The effect of small satel- 
lites, occasionally penetrating the disc at later times, can 
be seen in the third and fourth columns at L w r > 7 
and L w r rs 6 kpc, respectively. 

We note a strong variation of AL with cosmic time, 
where mergers dominate at earlier times (high z) and 
internal evolution takes over at t = 5 — 6 Gyr (corre- 
sponding to a look-back-time of ~6-7 Gyr, or z ~ 1). 
The latter is related to an increase with time in the bar's 
length and major-to-minor axes ratio as seen in the face- 
on plots, indicating the strengthening of this structure. 
Examining the top right panel of Fig. [TJ we find a con- 
tinuous increase in the bar's m = 2 Fourier amplitude 
with time, where the strongest growth occurs between 
t=4 and 8 Gyr. The effect of the bar can be found in the 
changes in angular momentum, AL, as the feature of neg- 
ative slope, centered on the CR (dotted-red vertical line), 
shifting from L s» 3.4 to «4.7 kpc. Due to the increase 
in the bar's amplitude, the changes in stellar guiding 
radii (vertical axis values) induced by its presence in the 
CR-region double in the time period 4.44 < t < 11.2 Gyr 
(bottom row of Fig. [T]). Until recently bars were not con- 
sidered effective at disc mixing once they were formed, 
due to their long-lived nature. We emphasize the im- 
portance of the bar in its persistent mixi ng of the inner 
disc th roughout the galactic evolution (see Minchcv et al. 
I2012allbl and discussion therein) . 

3. THE CHEMISTRY 

Our new approach is based on a detailed chemical evo- 
lution model for the thi n disc only (esse ntially the thin 
disc model presented bv iChiappinil [2009) . The idea be- 
hind this is to test if, once radial mixing is taken into 
account, we can explain the observations of both thin 
and thick discs without the need of invoking a discrete 
th ick-disc component (d i fferent from wha t was suggested 
in lChiappini et alJll997t IChiappinil 120091 ) . 

In the present model, the thin disc forms mainly via the 
slow accretion of extragalactic metal-poor gas (assumed 
to be of primordial composition). We assume the gas 
accumulates faster in the inner than in the outer regions 
(similar to what happens in the simulation discussed in 
Sec. [2]). Our code follows in detail a large number of 
chemical elements by properly taking into account the 
lifetime of stars of different masses (i.e., we do not use the 
insta ntaneous recycling approximation iChiappini et al.l 
Il997f h In the present work, we will concentrate only on 
Iron and OxvgerP^I. while other chemical elements will 
be discussed in the forthcoming papers of t his se r ies. 

As explained by IChiappini etaLl (|1997l . 120011 120031) , 
the main observational constraints of MW chemical evo- 
lution models are (i) the solar and present day abun- 
dances of more than 30 elements, (ii) the current stel- 
lar, gas and total mass densities at the solar vicinity 
(see Fig. [3J upper right panel), (iii) the present SFR 

10 We will often refer to the [a/Fe] ratios in the text, where in 
our case Oxygen is the chosen a element. 
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Fig. 2. — Properties of our detailed thin-disc chemical evolution 
model: Top, left: SFR in the solar neighborhood from our thin- 
disc chemical evolution model. Top, right: The total (solid), 
stellar (dashed) and gas (dotted) density evolution at the solar 
neig hborhood. The error b ars are the constraints from observations 
(see Chiappini et al. 2003, for references to the data). Bottom, 
left: SFR as a function of time; different curves correspond to 
different radii as indicated. Bottom, right: [Fe/H] gradients; 
different curves and colors correspond to different look-back times 
as indicated. 

(see Fig. [2 upper left panel), (iv) the present day su- 
pernovae rates of type II and la, and (v) the metal- 
licity distribution of G-dwarf stars (to be discussed in 
Sec. [5]). The SFH at the solar vicinity (see Fig. (2J upper 
left panel) is obtained by assuming a star formation law 
dependency on the gas density (Kennicut's law) and an 
exponentially decreasing infall with time, with an infall 
timescale r(R Si0 ) = 7 Gyr, where R„.^ = 8 kpc i s the 
Sun's galactocentric distance (see iChiappinU I2009L and 
references therein for details). From the chemical evo- 
lution constraints listed above, the one most affected by 
radial migration is the metallicity distributiorF 1 !. 

Note that this SFH is also in good agreement with the 
present amount of Deuterium measured in the ISM. The 
evolution of Deuterium has, apart from its obvious cos- 
mological interest, the role of a very sensitive marker of 
gas consumption in a give locale. All of the Deuterium 
in the Universe was created in the Big Bang and none 
has been released into the ISM ever since, but because 
it is consumed inside stars, its abundance steadily de- 
creases in direct proportion to the rate of star forma- 
tion. Thanks to the good knowledge of its primordial 
value from WMAP, a measure of Deuterium's present 
day abundance throughout the Galaxy provides an im- 
portant direct constraint on the stellar IMF, star forma- 
tion and infall histories at different locations (as these 
quantities define the levels of astration of Deuterium in 
the ISM). The th in disc model u s ed here is si milar to the 
one discussed by iRomano et al.l (|2003l |2006| ). and more 

Note that the age-metallicity relation, also affected by radial 
migration, was never used as a constraint in our chemical evolution 
models, due to the large scatter shown in the data. 



recently bv lLagarde et al.l (|2012l ). which well matches the 
pre-solar and present ISM measurements of D at the solar 
vicinity. 

Outside the solar vicinity, chemical evolution models 
are mainly constrained by the present day abundance 
gradients in the thin disc, and are thus more uncertain 
than the model for th e solar vicinity (see disc ussion in 
iChiappini etail I2001D . IChiappini etaLl (|200lD empha- 
sized that the radial abundance gradients predicted along 
the thin disc, as well as their time evolution, are strongly 
dependent on pre-enrichment and star formation thresh- 
old assumptions. Previous chemical evolution models in 
the literature considering the thin disc as an indepen- 
dent component could not anticipate such effects. Here 
we keep it simple by assuming neither threshold nor pre- 
enrichment; as a consequence, our models predict flat- 
tening of the gradients with time, similar to other pure 
thin-disc chemic al evolution models in the literature (e.g. 
IHou et aT]|2001|) . On the observational side, the question 
of whether the abundance gra dients steepen or flat ten 
with time is still not settled (see Stasinska et al. 20121 for 
a recent appraisal of the situation) . 

The chemical gradients we obtain are just a conse- 
quence of using different infall timescales at different 
radii. Here w e use the same r(r) expression as in 
IChiappini et al.l (|2001l ) down to r = 2 kpc. Chemical 
evolution models are often not computed for the inner- 
most 4 kpc to avoid dealing with the complete lack of 
observational constraints and with a region where differ- 
ent galactic components co-exist (bulge, bar, and inner 
disc). As here our main focus is the thin disc, we would 
like to assign chemistry down to r — and will make the 
conservative hypothesis that the chemistry computed for 
r = 2 kpc applies to all particles with r < 2 kpc, as well. 

The SFRs at different galactocentric radii are shown at 
the bottom left panel of Fig. [2j The resulting abundance 
gradients, for different look-back times, can be found in 
the bottom right panel of the same figure. 

Unfortunately, there are still large uncertainties also in 
the shape of present-day abundance gradients. From a 
recent compilation of t he be st data available in the liter- 
ature, [Stasihska^Tlil] (|2012t ) concluded that most claims 
of steepening or flattening of the present abundance gra- 
dients towards the centre or the outskirts of the MW are 
premature. The only result that seems reliable is the 
flattening of the disc planetary nebulae abundance gra- 
dient towards the centre of the Galaxy, as indicated by 
the comparison of the works by IChiappini et al.l (|2009l ) 
and I Henry et al.l (J2010). Whether this applies to other 
abun dance gradient trace rs is, as of yet, not established 
(e.g. iPedicelli et all 120091) . 

With our current SFH assumptions, the predicted 
present time gradient, shown by the black curve in the 
bottom left panel of Fig. [21 amounts to » —0.061 or 
w —0.057 dex/kpc, depending on the radial range we 
use for fitting (5-12 or 6-11 kpc, respectively). As we will 
see, slightly flatter values will be obtained (« —0.059 or 
w —0.058 dex/kpc) once radial migration is taken into 
account (Sec. It is beyond the scope of the present 
paper to check for the impact of other assumptions on 
the abundance gradients, but this will certainly be the 
focus of a future paper of this series. If the initial gra- 
dients turn out to be flatter than the ones adopted here, 
the effects of radial migration in some observables as the 
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age-metallicity relation will be, of course, weaker than 
the ones described in Sec. 15.41 

4. A NEW APPROACH FOR COMBINING CHEMISTRY 
AND DYNAMICS 

Despite the recent advances in the general field of 
galaxy formation and evolution, there are currently no 
self-consistent simulations, that have the level of chem- 
ical implementation required for making detailed pre- 
dictions for the number of ongoing and planned MW 
observational campaigns. Even in high-resolution N- 
body experiments, one particle represents ~ 10 4 ~ 5 solar 
masses. Hence, a number of approximations are neces- 
sary in order to compute the chemical enrichment self- 
consistcntly inside a simulation, leading to the several 
problems briefly discussed in the Introduction. Here, in- 
stead, we will assume that each particle is one staim 
and will implement the exact SFH and chemical enrich- 
ment from our chemical model (Sec. [3]) into our simulated 
galactic disc (Sec. [2]). This can be thought of as inserting 
the dynamics of our simulation into the chemical model, 
not the chemistry into the simulation. 

We start by dividing the disc into 300-pc radial bins. 
At each time output we randomly select newly born stars, 
by matching the SFH corresponding to our chemical evo- 
lution model at each radial bin. In case the simulation 
has formed more stars than needed, those are discarded. 
If we need more stars at a particular radius and time, we 
randomly select from the coldest stellar population to 
emulate the birth kinematics of the required number of 
extra stars. The particles selected at each time output 
are assigned the corresponding chemistry for that par- 
ticular radius and time. This is done recursively at each 
time-step of 37.5 Myr. 

In the just described manner, we follow the self- 
consistent disc evolution for 11.2 Gyr, selecting a tracer 
population possessing known SFH and chemistry enrich- 
ment. This bypasses all the problems encountered by 
previous chemodynamical models based on N-body sim- 
ulations. It also offers a way to easily test the impact 
of different chemical (or dynamical, see Appendix) pre- 
scriptions on the chemodynamical results. 

From a chemical point of view, there are two main 
simplifications in our approach. Firstly, we assume that 
radial flows do not play a substantial role in the global 
chemical enrichment of the Galaxy. This seems to be a 
good approximation as our differ ent chemical evolution 
zones are 2 kpc wide (see, e.g., ISpitoni and Matteuccil 
[20T1 . Secondly, we assume that stars do not contribute 
to the chemical evolution outside the zone where they are 
born, but either contribute only to the chemical enrich- 
ment within 2 kpc from their birth place, or never die. 
This assumption is valid for most of the stars because (i) 
the massive stars die essentially where they were born, 
due to their short lifetimes and (ii) low mass stars live 
longer than the age of the galaxy (never die) . We do not 
expect this simplification to affect our results by more 
than 10% for chemical elements made in low and inter- 
mediate mass stars, and even less for those coming from 
massive stars, as it is the case for Oxygen. 

From a dynamical point of view, the main simplifi- 

12 Dynamically, this is a good assumption, since the stellar dy- 
namics is collisionless. 



cations are (i) the resampling of the simulation SFH 
according to the chemical evolution model and (ii) the 
difference between the gas-to-total mass ratio expected 
from the chemical model and attained by the simulation. 
Although the discs in both the chemical and dynamical 
models grow inside-out, there are some offsets at particu- 
lar times. These differences in the SFHs are unavoidable, 
since the chemical and dynamical models are not tuned 
to reproduce the same star formation, although they are 
quite similar for most of the evolution (see also Sec. ??). 

As already mentioned in Sec. [U the disc gas-to-total 
mass ratio at the final time is ^0.12, consistent with 
the estimate of ~ 0.14 for the solar vicinity (Fig. [21 top 
right panel). However, at earlier times the gas fraction 
at the solar location, for example, can be overestimated 
by a factor of about 1.5-2 with respect to the chemi- 
cal model. One effect of a larger gas fraction can be 
a decreased suscep tibility to merger perturbations (e.g., 
iMoster et aLll2010fl . It may also give rise to a more un- 
stable disc, while suppressin g the bar instability (e.g., 
iBournaud and Combed [20021 ). Indeed, at earlier times 
we do find a weaker bar appearing more like a lens, 
which strengthens later, showing an increase in both 
its length and major-to-minor axes ratio (Fig. [T]). On 
the other hand, unreasonable disc instabilities due to a 
large gas fraction in the outer disc are not seen in our 
simulation. In fact, the spiral structure is about 15% 
of the background density and does not vary consider- 
ably, except when reinforced by mergers at earlier epochs. 
This amplitude is close to expectations for the MW (e.g, 
iDrimmel and Spergel|[200ll: iSiebert et al.ll2012t) . 

While the difference between the assumed (chemical 
model) and actual (dynamical model) stellar and gas 
densities can introduce some inconsistencies in the re- 
sulting dynamics, this would generally have the tendency 
of bringing more stars from the inner disc out, due to 
a larger bar expected at earlier times. As it will be- 
come clear later, a larger fraction of old stars coming 
from the bar's CR region to the solar vicinity would only 
strengthen our results. Overall, we do not anticipate the 
simplifications in our approach to affect significantly any 
of our results. In Sec. 19.11 we perform some tests which 
support this expectation. 

Overall, we do not anticipate the simplifications of our 
approach to affect significantly any of our results. On the 
other hand, what we gain with the above simplifications 
is a new tool to study the chemodynamics of our Galaxy, 
which is complementary to fully self-consistent models, 
and where the overall complex problem of galaxy assem- 
bly and evolution can be understood by pieces (i.e., same 
chemistry applied to different simulations and different 
chemistry applied to the same simulation). We antici- 
pate that this new approach will also be very useful to 
gain insights that can later be used in fully self-consistent 
simulations. 

5. THE RESULTING CHEMODYNAMICS 

After fusing the chemical and dynamical models as out- 
lined above, we are now in a position to (i) look for 
deviations from the predictions of the purely chemical 
evolution model and (ii) investigate the causes for these 
differences. 

To first illustrate the migration efficiency in our simu- 
lation, in the left panel of Fig. [3] we show the birth radii 
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Fig. 3. — Left: Birth radii of stars ending up in the "solar" radius (green bin) at the final simulation time. The solid black curve 
plots the total ro-distribution, while the color-coded curves show the distributions of stars in six different age groups, as indicated. The 
dotted-red and solid-blue vertical lines indicate the positions of the bar's CR and OLR at the final simulation time. A large fraction of 
old stars comes from the inner disc, including from inside the CR. Middle: [Fe/H] distributions for stars ending up in the green bin (left) 
binned by birth radii in six groups, as indicated. The total distribution is shown by the solid black curve. The importance of the bar's 
CR is seen in the large fraction of stars with 3 < ro < 5 kpc (blue line). Right: Density contours of the ro-[Fe/H] plane for local stars. 
The cyan curve shows our model solar-age metallicity gradient. Taking an error of ±1 dex in [Fe/H] we find a possible Sun birth radius of 
4.6 < ro < 7.6 kpc (where the horizontal green, transparent strip meets the cyan line). The imbedded histogram shows the density of stars 
in the likely ro-range, indicating a decline in the probability for ro < 5.5 kpc. 



of stars ending up in a solar neighborhood-like location 
(7 < r < 9 kpc, green bin) after 11.2 Gyr of evolution. 
The bar's CR and OLR are shown by the dotted-red and 
solid-blue vertical lines. The solid black line plots the 
total population, which peaks close to r = 6 kpc due to 
radial migration. An overdensity is seen just inwards of 
the bar's CR. The entire sample is also divided into six 
age-groups, shown by the curves of different colors and 
line-styles. The strongest effect from radial migration is 
found for the oldest stars (red curve), whose distribution 
has a maximum at r w 3 kpc, or inside the bar's CR. 
Note that locally born stars of all ages can be found in 
the solar neighborhood. A relatively smooth transition 
of the peak, from older to younger groups of stars, is 
observed; this is expected, since even for a constant mi- 
gration efficiency, older stars would be exposed longer to 
perturbations. While a wide range of birth radii is seen 
for all age groups, the majority of youngest stars were 
born at, or close to the solar neighborhood bin. 

5.1. The metallicity distribution 

The middle panel of Fig. [3] presents the metallicity dis- 
tribution of stars in our simulated solar neighborhood 
(green bin in left panel). In addition to the total dis- 
tribution (solid black line), six subsets of stars, grouped 
by birth radii, are shown by the different colors and line- 
styles. The locally born sample (green line) peaks slightly 
outside [Fc/H]=0, ends abruptly at [Fc/H]~ 0.15 and has 
an extended tail toward lower metallicities. This distri- 
bution is the same as the one predicted by our pure-thin- 
disc chemical evolution model with the SFH shown in the 
top left panel of Fig. [2] Our chemo-dynamical model in- 
stead, predicts that the majority of migrators comes from 
the inner region (3 < ro < 5 and 5 < ro < 7 kpc bins), as 
expected. A smaller, but sizable, fraction of stars orig- 
inating outside the solar neighborhood is also observed 
(orange and red curves). 

While the low-metallicity end in the total distribution 



comes from contribution by all initial radii, the tail at 
larger metallicities extending to [Fe/H]~ 0.6 results ex- 
clusively from stars with 3 < ro < 5 kpc (note that the 
5 < ro < 7 kpc bin contributes up to [Fe/H]~ 0.25 only). 
This is the region just inside the bar's CR, which is where 
the strongest outward radial migration occurs, as we dis- 
cussed in Sec. [5]and showed in the bottom row of Fig. [TJ 
At s maller radii stars are trapped in the X\ b ar orbits 
(e.g.. IContopoulos and P apavan n"opouloslll980[ ). making 
outward angular momentum transfer difficult. Neverthe- 
less, stars originating from 1 < ro < 3 kpc still appear 
in the local bin, with a fraction comparable to that of 
objects coming from ro > 11 kpc. We discuss how the 
metallicity distribution relates to observations in Sec. 16.21 
(Fig. HO}. 



5.2. The birth place of the Sun 

With the care taken into defining a proper solar radius 
in a simulation with MW characteristics, combined with 
MW chemistry, it is tempting to make an estimate for 
the solar birth radius, which emerges naturally from our 
chemodynamical model. 

Already from the left panel of Fig. [3] we can see that 
the age distribution to which the Sun belongs (dotted- 
cyan curve), suggests a most likely solar birth location at 
r « 6 kpc and similar probability of it being born in situ 
or at r « 5 kpc. By including the chemical information 
we can further improve this estimate. 

The right panel of Fig. [3] displays density contours 
of the ro-[Fc/H] plane for all local stars. The cyan 
curve shows our input solar-age (4.6 Gyr look-back time) 
metallicity gradient. Assuming an error of ±1 dex in 
[Fe/H], we find a possible Sun birth radius in the range 
4.4 < ro < 7.7 kpc (where the horizontal green, trans- 
parent strip meets the cyan curve). 

The imbedded ro-histogram in the right panel of Fig. [3J 
shows the density of stars in the likely ro-range, esti- 
mated from stars in a narrow age-bin around the solar 
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Fig. 4. — Top: The left panel plots age versus [Fe/H] for different radii, resulting from our input chemical model (left). The middle panel 
shows stellar density contours of the resulting relation after fusing with dynamics, for the "solar" radius (7 < r < 9 kpc). The overlaid 
lines show the input chemistry for the same radii as in the left. The right panel plots the metallicity distributions for different age bins. 
Bottom: Same as above but for the [Fe/H]-[0/Fe] relation. There is some contribution from stars born ~ 2 kpc and hardly any from 
r > 14 kpc, consistent with the birth radii distributions shown in Fig. \3\ 

due to the predominance of migrators coming from the 
inner disc, as seen in Fig. [3J Similarly, contours below 
the cyan curve result from stars coming from the outer 
discQ The bump at age~ 9.5 Gyr and [Fe/H]« -0.1 re- 
sults from the strong, merger related migration efficiency 
at that time (see Sec. 16. 3p . Identifying such a structure 
in the observed age-metallicity relation would be indica- 
tive of a strong change in the migration efficiency with 
time (due to the 1:5 mass ratio merger ~ 9 Gyr ago in 
this case). 

The right top panel of Fig. 2] shows the metallicity dis- 
tributions for different age bins. This is consiste nt with 
recent observational results fe.g.. lHavwoodll2012l and ref- 
erences therein) , in showing that the scatter in metallic- 
ity increases towards older ages. We find the following 
dispersions in [Fe/H] for age groups from young to old: 
0.11, 0.14, 0.16, 0.17, 0.15, and 0.7. In all but the oldest 
bin the distributions appear gaussian. The shape of the 
red histogram (oldest stars) can be approximated better 
by a log-normal distribution because of its strong nega- 
tive skewness. The metallicity peak of this oldest popu- 
lation is around —0.5, with a r ange of values in very g ood 
agreement with the results bv lSoubiran et all (|2003f ) for 
their thick-disc selection of clump giants. As we found 
our from Fig. [3J more than 50% of these old stars were 

13 The deficiency of stars seen at age> 9 Gyr ([Fe/H]Ri —0.8) in 
the lowest density contours is an artifact due to the slightly lower 
SFR in our simulation at r > 10 kpc; note that this only affects a 
negligibly small number of the oldest stars at this large radius. 



201 If ). We find the highest probability to be around 
5.6 kpc, followed by 7 kpc. Note that this estimate is 
dependent on the migration efficiency in our simulation 
and the adopted chemical evolution model. 

5.3. The age-[Fe/H] and [Fe/Hj-[0/FeJ relations 

We showed in Figs. [T] (bottom) and [3] (left) that radial 
mixing is significant in the simulation we use here. It 
is therefore, very interesting to find out how much the 
age-metallicity (AMR) and [Fe/H]-[0/Fe] relations are 
affected. 

To understand better our results, in the first row, left 
panel of Fig. 2] we show stellar age versus [Fe/H] coming 
from our input chemical model, for different radii, as in- 
dicated in the lower panel. The middle panel shows stel- 
lar density contours of the resulting relation after fusing 
with dynamics, for the "solar" radius (7 < r < 9 kpc). 
The overlaid curves highlight the input chemistry for the 
same radii as in the left, giving an insight about the ori- 
gin of stars found in this localized bin. We see that stars 
with ro > 14 kpc do not enter the solar bin. How much 
the AMR of locally born stars flattens due to radial mi- 
gration can be seen by comparing the cyan curve, which 
corresponds to 8 kpc, to the highest density contours 
(blue and black). Some flattening is observed, mostly 
for ages >5 Gyr, however the final distribution is by no 
means flat. The excess of stars at high metallicities is 
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born at r < 5 kpc. The causes for their outwards migra- 
tion is discussed in Sec. 16.31 

The bottom row of Fig. [4] is the same as the top one, 
but for the [Fe/H]-[0/Fe] relation. The relationship be- 
tween [O/Fe] and age, displayed in the rightmost panel, 
shows good correlation between the two, assuring that 
a-elements can be a good age-indicator. However, this 
figure also shows that the correlation is best seen for older 
ages. For age<6 Gyr, the different distributions overlap 
significantly (see also Fig. [51 bottom). This should be 
kept in mind when drawing conclusions from large spec- 
troscopic surveys, where often the [a/Fe] ratio is used as 
a proxy for age. 

By correcting for the spectroscopic samp ling of stel- 
lar sub -populations in the SEGUE survey, IBovv et al.l 
(|2012bf ) showed recently that the bimodality seen in the 
[0/Fe]-[Fe/H] distribution for the uncorrected sample 
disappears. As the middle bottom panel in Fig.Ulshows, 
our results are in agreement in that our simulated, unbi- 
ased [0/Fe]-[Fe/H] stellar density distribution is smooth. 
In Sec. 16.11 we show that biases introduced by the pref- 
erential selection of thin- and thick-disc stars, as com- 
monly employed in observations, can result in a gap in 
the [0/Fe]-[Fe/H] plane. 

5.4. Migration effects on the chemical gradients 

In Fig. [5] we show how the initial chemical gradients 
are affected by the disc dynamics. We divide the galac- 
tic disc into six age-bins and plot the [Fe/H] and [O/Fe] 
gradients by lines of different colors, as indicated in the 
bottom panel. Solid and dotted line-styles represent the 
initial (chemical) and final (chemodynamical) states, re- 
spectively. The bar's CR and OLR (estimated at the 
end of the simulation) are given by the dotted-red and 
solid-blue vertical lines, respectively. A strong flattening 
in [Fe/H] is seen for the older populations, however, the 
younger stars are much less affected. Remarkably, for 
stars younger than 2 Gyr the final gradient is very simi- 
lar to the initial one out to ~ 12 kpc. This is related to 
the fact that the majority of stars in this age-bin are born 
near their final radii, as we found out was the case for the 
"solar" radius (see dotted black line in Fig. [3l left). Note 
that the bar's CR acts as a pivot point around which the 
initial metallicity profiles turns counterclockwise, raising 
values in the disc and lowering them in the bulge/bar 
region. This is easy to understand by recalling that the 
bar CR radius is where the strongest exchange of angular 
momentum occurs (see Fig. [TJ bottom row) and is one 
example of the importance of considering the effect of 
the bar in modeling the Galactic disc. 

To give a quantitative idea for the change in the metal- 
licity gradient of young stars, we fitted lines to the solid 
(initial) and dotted (final) black lines in Fig.[5j top panel. 
We estimated that the original gradient of « —0.061 or 
?» —0.057 dex/kpc, depending on the radial range we 
use for fitting (5-12 or 6-11 kpc, respectively) changes to 
w —0.059 or s» —0.058 dex/kpc at the end of the simu- 
lation. This small effect on the gradient of the youngest 
stellar population is reassuring, since this was used as a 
constraint for our input chemical evolution model. 

Incidentally, it turns out that the CR (which is consis- 
tent with what is expected for the MW) is almost exactly 
at the radius where the input metallicity gradient begins 
to steepen (Fig. 03 top). While the chemical model we use 
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Fig. 5.— The effect on the initial [Fe/H] (top) and [O/Fe] (bot- 
tom) gradients for different stellar age groups. The solid and dot- 
ted color curves show the initial and final states, respectively. Note 
that, while strong flattening is observed for the older populations, 
the metallicity gradient for the youngest stars {age < 2) is hardly 
affected at r < 12 kpc, thus justifying the use of our chemical 
model, which uses this as a constraint. 

was not motivated by the effect of the bar and is, in fact, 
a somewhat extreme case regarding the strong increase 
in [Fe/H] at r < 5 kpc, this choice may well be the cor- 
rect one, considering the stronger chemical enrichment 
expected inside bars (e.g.. lEllison et al1l2011h . 

Much weaker flattening is seen for the [O/Fe] profiles, 
related to the progressively weaker radial variation of 
older samples. Due to the reversed gradients compared 
to [Fe/H], a clockwise turn is now seen around the CR 
"pivot point" . In both panels of Fig. [5] the values in 
the innermost disc region are unaffected. This is be- 
cause most of these stars are trapped inside the bar, as 
mentioned above, preventing their escape once the bar is 
formed. 

The prediction that the [O /Fe] radial profiles are essen- 
tially preserved is an important result and can provide a 
powerful constraint on different input chemical evolution 
models. The shapes of the [0/Fe]-profiles for different 
ages bins are a direct consequence of our adopted varia- 
tion of the infall-law with radius (and hence the SFHs at 
the different positions). In the near future, these would 
be possible to measure by combining the g ood distances 
and a ges expected from the CoRoT mission (Bag lin et al.l 
2006), with abundance ratios obtained by spectroscopi c 
follow-up surveys (see iChiappml [20121: iFreemanl feO^). 
For the young populations, this should be already possi- 
ble to be obtained from the observations of open clusters, 
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Fig. 6. — Left column: The age-velocity relations (AVR) result- 
ing from our model in the radial (top panel) and vertical (middle 
panel) directions. The black solid lines show the AVR estimated 
from all stars in the annulus 7 < r < 9 kpc. The green, red and 
blue squares plot the locally born stars, stars born in the outer 
disc (rrj > 9 kpc), and in the inner disc (ro < 7 kpc), respectively. 
No significant contribution is found from samples arriving from 
different galactic radii. The bottom panel shows the relationship 
between age and [O/Fe] for different birth-radius bins. Right col- 
umn: The radial and vertical velocity dispersions as functions of 
[O/Fe]. The reason for the flattering at [O/Fe]<-0.1 is the mixture 
of stars of age< 4 Gyr. This can be seen easily in the bottom panel, 
which is similar to the one on the left, but with the horizontal and 
vertical axes reversed. The orange lines in the top row show he 
stellar rotational velocity, V^. 

e.g., with t he ongoing; Gaia-ESO (jGilmore et al.ll2012T) or 
APOGEE (jMaiewski et aLllMTol ) surveys. 

5.5. The age-Velocity and [O/Fe] -Velocity relations 

An important constraint in MW disc modeling is the 
observed velocity dispersion in the solar neighborhood. 
Due to the lack of good age estimates, the shape of the 
age-velocity relation (AVR) has been a matter of de- 
bate. While a powe r-law of the form a = t°- 2 -°- 5 has 
been proposed (e.g . , iWielenl Il977t iDehnen and Binnevl 
\199§. iBinnev et alJ I2000D . some works have predicted 
an increase for stars younger than 2 Gyr, followed by 
a saturation up to ~10 Gyr, followed by a strong in- 
crease associated with the thick disc, possibly related 
to strong merger acti vity (or a single major merger) 
at high redshift (e.g.. [Freemanl Il991t lEdvardsson et al.l 
H993 iQuillen and GarnettJl200oh . 

In the left column, top two panels of Fig. [5] we present 
the radial and vertical velocity dispersions, ay and a z , 
as functions of age for stars in the simulated solar neigh- 
borhood (7 < r < 9 kpc). The solid black lines plot the 
total sample, while the differently color-coded squares 
show the locally born population (green), stars born in 
the outer disc (red) and stars born in the inner disc 



(blue). Examining the total population, we note that 
for ages greater than ~ 3 Gyr, oy flattens significantly, 
rises only by 5 km/s in the time interval 3-8 Gyr, after 
which the slope increases again, finally reaching 60 km/s 
for the oldest stars. The a z behavior is somewhat differ- 
ent, with a smooth growth from to 8 Gyr, followed by 
a steeper increase, even more drastic than ay, reaching 
about 50 km/s for the oldest stars. This is roughly con- 
sistent with the second scenario mentioned above. We 
must note, however, that some increase in velocity dis- 
persion with age is always seen, except possibly in the 
range 6-8 Gyr for a r . 

We now compare the contribution to the AVR from 
stars born in situ, inside the "solar" circle, and out- 
side it. These three subpopulations, indicated by squares 
of different colors, show remarkably similar AVRs. The 
subsamples born outside/inside the radial bin considered 
have slightly lower/higher velocity dispersion as a func- 
tion of age (except at the oldest ages, se e below). This 
is consistent w ith the recent findings by Minchev et al.l 
(|2011bl I2012H) . where it was demonstrated that migra- 
tors are generally not expected to heat the disc, due to 
the approximate conservation of the average radial and 
vertical actions in the absence of external perturbers. It 
is remarkable that we do not see any significant contribu- 
tion here either, although the disc we study is perturbed 
by satellites, thus breaking the action invariance. We 
note that no massive mergers have occurred in our sim- 
ulation for the last ^8 — 9 Gyr, which could explain 
this behavior. Interestingly, we find that, for the old- 
est population, the locally born stars have substantially 
higher velocity dispersions (especially a z ) compared to 
the subsample coming from the inner disc (contrast green 
and blue squares for age>10 Gyr in the second row, left 
panel). This is in drastic contrast to the assumption 
made in a number of studies, that stars coming from the 
inner disc should heat vertically. The reason that, for the 
oldest population, the hottest stars were born in situ or 
outside the "solar" circle is related to the stronger effect 
of merger-induced perturbations in the outer disc (due 
to its exponential density decrea se with radius), gener - 
ally resulting in disc flaring (e.g.. IBournaud et al.l l2009). 
The larger ovvalues for stars of age> 10 Gyr coming 
from the outer disc, compared to the locally born pop- 
ulation, supports this conclusion (see also Sections 16.31 
and US). 

Due to the lack of good age estimates, the [a- 
elements/Fe] ratios, such as [O/Fe], are often used as 
proxies of age. To see how well age- velocity and [O/Fe]- 
velocity relations match, in the top two panels of the 
right column of Fig. [6] we plot the radial and vertical ve- 
locity dispersions as functions of [O/Fe]. In agreement 
with the bottom right panel of Fig. [J [O/Fe] appears to 
follow age closely at the a-rich end, but not that well at 
low values. The stars born locally and outside the bin 
considered (green and red squares), show no information 
for [0/Fe]< -0.07 and [0/Fe]< 0, respectively. This 
can be related to our chemical evolution model (see the 
bottom left panel of Fig. SJ) . 

What gives rise to the flattening at [0/Fc]< —0.1 seen 
for both ay and a z ? By inspecting the bottom panel of 
Fig. 03 we see that stars with [0/Fe]< —0.1 are born at 
tq < 7 kpc and have age< 8 Gyr, correlating inversely 
with birth radius. Therefore, the majority of stars with 
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[0/Fe]< —0.1 currently in the local bin, should have 
ages< 4 Gyr (see also bottom panel). Such a mixture of 
ages would naturally result in higher velocity dispersions. 
Indeed, the mean ay and a z for stars of age< 4 Gyr seen 
in the left column of Fig. [SJ match well the correspond- 
ing, approximately constant, values at [0/Fe]< —0.1 in 
the right column. 

Concentrating on the local population (green squares) , 
we note markedly lower velocity dispersion values, com- 
pared to the total sample, at the lowest [O/Fe] achieved. 
Inspecting the bottom panel of Fig. [5l it is evident that 
only the youngest locally born stars possess such low 
[O/Fe], thus the velocity dispersions are similar to those 
of the youngest stars seen in the left panel. A similar ar- 
gument explains the low oy and u z -values at [0/Fe]< 0.1 
for stars originating in the outer disc (red squares). As 
in the plots on the left, for the high-[a/Fe] stars we find 
higher velocity dispersion values for locally born stars 
and inward migrators, compared to stars born in the in- 
ner disc. The reason for this was outlined above. 

The orange curves in the top row of Fig. [6] show the 
stellar rotational velocity, V^, as a function of age (left) 
and [O/Fe] (right). An inverse correlation between V$ 
and ay is expected, due to the asymmetric drift effect. 
Naturally, a constant decrease in is seen with in- 
creasing age. Note that the oldest stars lag the younger 
population by ~ 50 km/s, similar to the rotational lag 
found for the t hick-d isc selection of clump giants by 
ISoubiran et al.l (|2003[ ). for example. The more compli- 
cated behavior found in the right panel is expected in 
view of our discussion above, where at low [O/Fe] the 
population is a mixture of stars with ages up to ~ 4 Gyr. 
An upcoming paper is dedicated to the interpretation of 
the metallicity-V^ relation. 

5.6. The relation between age and [O/Fe] 

We now briefly examine the relation between age and 
[O/Fe] for stars ending up in our simulated solar vicin- 
ity. In the bottom left panel of Fig. [6] we plot [O/Fe] 
as a function of age for six birth-radius bins at the final 
simulations time. The pink line shows the total sample, 
which, for age> 4 Gyr, is seen to follow the stars with 
birth radii 5 < ro < 7 kpc (cyan) , rather than the locally 
born population (green). This goes well with Fig. [31 left 
panel, where it was seen that the ro-distribution peaks 
at around 6 kpc. For age< 9 Gyr, approximately lin- 
ear relations are found for stars born in the same radial 
bins. However, due to contamination from radial migra- 
tion, for the total sample the scatter increases strongly 
with decreasing age, in agreement with the bottom right 
panel of Fig. |4j 

The bottom right panel of Fig. [6] shows age as a func- 
tion of [O/Fe]. It is now easy to see the contribution 
from different galactic radii to the AVRs in the panels 
above. As we found out earlier by different reasoning, 
stars at [0/Fe]< —0.1 arc a mixture of ages< 4 Gyr. 
While following stars born at ro ~ 6 — 8 kpc (green and 
cyan bins) at age> 3.5 Gyr, the total sample (pink line) 
shows an abrupt change for age< 3 Gyr, shifting toward 
the bin of birth radii 3 < ro < 5 kpc. This introduces a 
strong non-linearity in the age- [O/Fe] relation for these 
younger ages. 

In summary, the approximately linear relation between 
age and [O/Fe] for stars born at the same radius is de- 



stroyed when radial migration is taken into account. 
When grouped by common birth radii, each sample is 
a quite good age-indicator, in the covered [O/Fe] region, 
e.g., in the range —0.1 <[0/Fe]< 0.4 for the locally born 
stars. However, due to the inevitable radial mixing, the 
overlap of these sub-populations, combined with the off- 
set in the [0/Fe]-range covered, contaminate this rela- 
tion. 

5.7. The vertical structure of the local disc 

We now study the density and vertical velocity as 
functions of height above the disc plane in the annulus 
7 < r < 9 kpc. Unlike the common division into thin 
and thick discs, we chose to divide our sample in groups 
by age, [O/Fe], or [Fe/H]. 

In the first row of Fig. [7] we show the vertical velocity 
dispersion as a function of height above the disc plane, 
a z (z), for stars in the "solar" vicinity (7 < r < 9 kpc) 
binned in six different groups by age (triangles, left), 
[O/Fe] (squares, middle), and [Fe/H] (squares, right). 
The corresponding bin values are shown to the right 
of the figure. Very similar increasing trends are seen 
for young stars, low-[a/Fe] and metal-rich bins [a z ~ 
10 — 20 kpc). Interestingly, as populations get older, 
relations flatten around ~ 50 km/s for the oldest /high- 
[a/Fe] /metal-poor samples. Both these trends and val- 
ues are similar to the recent observational results by 
IBovv et all l|2012d) . We note that in our work sub-pop- 
ulations are in narrow bins of [O/F e] or [Fe/H] , but no t 
both at the same time, as done by IBovv et al.l (|2012cf) . 
Nevertheless, the similarity with their Fig. 2, bottom, is 
obvious. 

The second row of Fig.[7]displays the vertical disc scale- 
height at the "solar" radius binned in the same groups 
as in the first row. We find that single exponentials can 
fit reasonably well each subsample, as well as the total 
population. The color-coded values, h, indicate single 
exponential fits, over-plotted by the dotted lines (dashed 
black line for the total population). No smoothing is 
applied. Gradual trends are seen for all panels, where 
values of h grow with increasing age, increasing [O/Fe], 
and decreasing [Fe/H]. 

Older samples are always expected to have higher ve- 
locity dispersion, on the average, than younger ones, due 
to their being exposed longer to a variety of heating 
agents. As older samples are high-[a/Fe] and metal-poor, 
the trends seen with mean height above the disc plane are 
not surprising. It is not necessarily obvious why single 
exponential would provide good fits to stars in narrow 
abundance or age bins. As we show in Paper II, this is 
not expected throughout the entire disc extent. 

The variation of scale-height for different sub-popu- 
lations, as well as the fact that we can fit single ex- 
po nejrtials_to them, are consistent with the findings 
of IBovv et all (l2012ah . based on G-type SEGUE dwarf 
data. Again, we point out that we have constrained sam- 
ples either by [O/Fe] or [Fe/H], but not both, as Bovy et 
al. did. 

It is also interesting to consider the variation of [Fc/H] 
with height above the disc plane, z, as well as with mean 
vertical velocity, v z . In Fig. [S] we show these trends for 
the same bins of age (left) and [O/Fe] (right) used for 
Fig. [7] The solid black lines in each panel display the 
total number of stars in the annulus 7 < r < 9 kpc. We 
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Fig. 7. — First row: The vertical velocity dispersion as a function of height above the disc plane, cr z {z), for stars in the "solar" vicinity 
binned in six different groups by age (left), [O/Fe] (middle), and [Fe/H] (right). The corresponding bin values are shown to the right of the 
figure. Very similar increasing trends are seen for young stars, o-poor and metal-rich bins (<r z ~ 10 — 20 kpc). As populations get older, 
relations flatten. Second row: Same as above, but for the vertical scale-heights. The color-coded values, h, indicate single exponential 
fits, shown by the dotted lines (dashed black line for the total population). Smooth increase of scale- height is seen with increasing age, 
increasing [O/Fe], and decreasing [Fe/H]. 



consider stars above and below the plane and take the 
mean absolute values \z\ and \v z \. A negative trend is 
seen for both \z\ and \ v 7 \ as expected from ob servations 
(e.g., iKatz et al.ll20TH iSchlesinger eTaT]l20lih . 

When divided into age or [O/Fe] subsamples, the vari- 
ation of [Fe/H] is mostly lost for each individual group. 
This is true for both \z\ and \v z \. It is then clear, that 
the downtrend in the total sample comes from the fact 
that younger (or a-young) stars are, on the average, both 
metal-rich and more numerous at small distances from 
the plane. With increasing height the density of older 
(large [O/Fe] -ratios), and thus more metal-poor, pop- 
ulations continuously increases, resulting in the overall 
negative gradient. 

Similarly, younger stars are confined close to the disc 
and, thus, have relatively low vertical velocities. With in- 
creasing age (and thus velocity dispersion) the amplitude 
of vertical oscillations increases, while the stellar density 
decreases. This, again, results in the negative trend of 
[Fe/H]. The flattening at \v z \ > 90 km/s is due to the 
fact that the number of stars with such high velocities is 
almost exclusively dominated by the oldest/ [0/Fe]-rich 
stars. 

The approximately flat mean \z\ — [Fe/H] and 
|ujt|— [Fe/H] relations for stars of narrow age- or 
[0/Fe]-bins is an extension to the similar behavior seen 
in the \z\ — a z ^ r plots above and provides additional 
predictions for Galactic surveys. 

6. RELATION TO RECENT OBSERVATIONAL RESULTS 

We have already shown that a number of the basic 
observed chemistry and kinematics of solar neighborhood 
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Fig. 8. — First row: [Fe/H] as a function of height above the 
disc plane (absolute value), \z\, for stars in the "solar" vicinity. 
Different colors indicate age (left) or [O/Fe] (right) for the same 
bins as in Fig. [4] Second row: Same as above but the vertical 
axes show the absolute mean vertical velocity, \v z \. The solid black 
curve in each panel shows the relation for the total sample. 



stars can be explained by our model. We now perform 
some further tests, accounting roughly for sample biases 
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Fig. 9. — Selection effects can result in a bimodality in the 
[Fe/H]-[0/Fe] plane. The top panel show the unbiased stellar den- 
sity distribution, as in Fig. [4] In the bottom plot we have applied 
the selection criteria used by Bcnsby ct al. (2003). 

encountered in different surveys. 

6.1. Is the bimodality in the [Fe/H]-[0/Fe] plane due to 
selection effects? 

Starting with the assumption of two distinct entities - 
the thin and thick discs - many surveys s elect stars ac- 
cording to certain kinema tical criteria (e.g., 
20031 IReddv et all [20061) . For example, 
( 2003) defined a (now widely used) method for prefer- 
entially selecting thin- and thick-disc stars with a prob- 
ability function purely based on kinematics. 

Starting with the smooth [Fe/H]-[0/Fe] stellar density 
distribution, presented in Fig. [H we now employ the same 
technique as lBensbv et al.l (|2003[ ) and extract a thin- and 
thick-disc selections from our initially unbiased sample. 
In order to obtain a similar number of thin- and thick- 
disc stars, as done in most surveys, we randomly down- 
sample the thin-disc selection. 

The results are presented in Fig. HI where the top panel 
shows the unbiased distribution. In the bottom panel we 
have applied the selection criteria as described above. We 
find a discontinuity at [O/Fe]«0.15-0.2, similarly to what 
is seen in a number of observational studies. This test 
shows that selecting stars kinematically can give prefer- 
ence to the two extremes of the distribution (oldest vs 
youngest), resulting in the gap we find. 

6.2. Changes in the MDF for samples at different 
distances from the disc plane 

To be able to compare our results to different MW 
surveys we need to constrain spatially our sample. For 
example the GCS stellar sample is confined to within 
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Fig. 10. — Changes in the metallicity distribution function 
(MDF) for samples at different distances from the disc plane. First 
column: Distribution of birth radii for stars ending up in the bin 
7 < r < 9 kpc (green rectangle) at the final time. The bar CR 
and OLR are indicated by the dotted red and solid blue vertical 
lines. From top to bottom, we have applied the selection criteria 
\z\ < 0.5, \z\ > 0.5, and \z\ > 1.5 kpc. Second column: MDF 
for the same particle subsets as on the left. The black line show 
the total population for each selection, while the red and blue col- 
ors indicate [0/Fe]> 0.3 and [0/Fe]< 0.3, respectively. A shift in 
the peak from [Fe/H]=0 to negative values is seen with increasing 
distance from the disc plane. 



100 pc from the Sun, while the SEGUE and RAVE sam- 
ples explore much larger regions, but miss the local stars. 
To see the biases introduced in the MDF when distance 
cuts are applied, we explore three possibilities. 

Figure [TU] shows our predicted metallicity distribution 
function at the solar vicinity, for different heights from 
the plane. From top to bottom, we have applied the se- 
lection criteria \z\ < 0.5 (which is similar to what is seen 
by GCS), \z\ > 0.5 , and \z\ > 1.5 kpc (sim ilar to what is 
seen in SEGUE bv lSchlesinger et al.ll20TTI) . The left col- 
umn shows the distribution of birth radii for stars ending 
up in the bin 7 < r < 9 kpc (green rectangle) at the fi- 
nal time, whereas the right column shows the metallicity 
distribution for the same particle subsets as on the left. 
In the right column the black line show the total pop- 
ulation for each selection, while the red and blue colors 
indicate [0/Fe]> 0.3 and [0/Fe]< 0.3, respectively. It 
can be seen that the MDF peak shifts from solar at low 
heig hts from the plane (cons istent with the GCS sample, 
e.g.. ICasagran de et al.|[20Tll) . to a value below solar, as 
we select only particles with heights above 500 pc (con- 
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Fig. 11. — First row: [Fc/H]-[0/Fe] stellar density distributions 
for metal-poor, high-[0/Fc] stars ending up in the cylinder 7 < r < 
9 kpc, 0.5 < \z\ < 2.5 kpc. Stars are separated into a "cold" and 
"hot" samples with corresponding velocity dispersions a r 65 
and o> ss 28 km/s. Second row: Birth radii distributions for 
two age groups, as indicated. Third row: Same as the second 
row, but at t = 3 Gyr, i.e., 0.5 Gyr after the most massive merger 
has concluded. A large fraction of stars for both the hot and cold 
samples are in place at that time. 

sistent wit h RAVE - Boeche et a l. 2012, In prep, and 
SEGUE - ISchlesinger et alll20ll . In the latter case a 
value of [Fe/ H]~ —0.3 is found, strikin gly similar to the 
one found by ISchlesinger et al.l (|2011D . Note that even 
at large distance from the plane, the MDF is dominated 
mostly by stars with [0/Fe]< 0.3, whereas the higher 
[O/Fe] values are always contributing to its metal-poor 
tail. This sugge sts a bias in the d istributions obtaine d by 
iLee et al.1 (|2011[ ). as discussed bv lBovv et all (|2012al lbl|cl). 

6.3. The origin of stars with low-metallicity and 
high- [O/Fe] ratios 

Analyzing the SE GUE G-dwarf data, 
iLiu and van de Venl ([2012) argued that (internally 
driven) radial migration cannot be responsible for 
the origin of the hot, metal-poor, high- [O/Fe] stars 
and proposed that they have formed through early-on 
gas-rich mergers. 

We now search for the origin of the metal-poor high- 
[O/Fe] stars in our simulation. To approximate the spa- 
tial coverage of the SEGUE G-dwarf sample, we consider 
stars in the annulus 7 < r < 9 kpc and distance from the 



plane 0.5 < \z\ < 2.5 kpc. As in ILiu and van de Venl 
i|2012ft (their group III) we consider stars with —1.2 < 
[Fe/H] < -0.6 and 0.25 < [O/Fe] < 0.47. Finally, we 
apply a cut in eccentricities at e = 0.2, thus separat- 
ing our sample into a "cold'Q and "hot" subgroups with 
corresponding radial velocity dispersions ay ~ 27 and 
67 km/s. The ratio of cold to hot orbits is ~ 0.3. 

In the top row of Fig. [11] we plot stellar density con- 
tours of the [Fe/H]-[0/Fe] relation for the hot (left) and 
cold (right) orbits. As in Fig. |4j we overlay the input 
chemistry for different initial radii, where the thick cyan 
curve corresponds to r — 8 kpc and the rest are sepa- 
rated by two kpc (black for r = 2 kpc). We can already 
see that a large number of stars come from the inner disc, 
especially for the hot sample. 

The second row shows that initial radii span the range 
1 < ro < 10 kpc for both samples, but the distribution 
for the cold orbits is heavier near the solar vicinity. The 
majority of stars have age> 10 Gyr (solid red line) with 
initial radii spread all over the disc. A much smaller 
group of stars with 9 < age < 10 Gyr (dotted orange) 
originates at the solar vicinity. The ro-distribution of 
stars older than ~10 Gyr is an example of extremely effi- 
cient migration, usually not seen as the effect of internal 
perturbations. 

To get an idea of when these stars arrive to the solar 
neighborhood, in the third row of Fig. [TT] we plot the 
distribution of initial radii obtained 3 Gyr after the be- 
ginning of disc formation. At this early stage we already 
see distributions similar to the ones we find "today" . For 
example, there are three peaks in the ro-distribution at 
?» 3.5, 5, and 7 kpc for both hot samples and the cold 
selection today. Note that not all of the stars found in 
the solar vicinity at t = 3 Gyr would remain there for 
the next ~ 8 Gyr due to the continuous migration. Ad- 
ditionally, part of the cold sample at t = 3 Gyr would 
heat up, ending up in the hot one we see today. 

What causes the extreme migration at t < 3 Gyr? As 
mentioned in Sec.[2l the last massive merger (1:5 mass ra- 
tio) begins at t » 1.5 and is disrupted by t « 2.5 Gyr, or 
just before the t — 3 Gyr distributions shown in Fig. [TTJ 
The strong effect on the changes of stellar angular mo- 
menta around this time were presented in Fig. [TJ bot- 
tom left panel, leaving no doubt that this event is re- 
sponsible for the vehement migration we find. In the 
Appendix (Sec. 19. 3| Fig. [T4|) we demonstrate that the 
merger-induced migration suppresses heating both in the 
vertical and radial directions. On the contrary, after the 
merger has concluded, the already preheated inner disc 
is a source of hot stars, thus contributing to disc heat- 
ing via migration induced by the bar and spirals. As the 
disc grows further from gas accretion, this effect becomes 
much smaller. 

Further contribution from the bar in migrating old 
stars out is expected at t > 3 Gyr, as inferred from the 
increase in the ro-distributions at ~ 3.5 kpc at the fi- 
nal time (middle row). It should be kept in mind that 
migration effici ency decreases with increasin g velocity 
dispersion (e.g.. lComparetta and Qui llcn 2012), present- 
ing an extreme case in the current kinematic selection. 

14 Here we are somewhat stretching the meaning of "cold" orbits. 
Usually, for stars at the solar radius, this would be a sample with 
o r ~ 5 — 10 km/s. 
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Therefore, the bar and spirals at later times would be 
mostly transferring cold orbits, populating the region at 
r q 5, 6.5 kpc in the co l d distr ibution. 

iLiu and van de Venl (|2012t ) argued that the old, hot, 
metal-poor, high-[0/Fe] stars were born hot in gas- rich 
mergers. Our analysis above suggests a very similar in- 
terpretation (although our disc does not experience a 
large number of gas-rich mergers), in the sense that stars 
were born hot and further heated by mergers, which also 
gives rise to large-scale, very efficient migration. We 
also find, however, that further migration at later times, 
driven by the bar and spirals, is essential, as discussed 
above. 

To see when the oldest stars migrate to the solar neigh- 
borhood in the lack of early massive mergers, we per- 
formed a different model realization for which the dy- 
namics is quiescent and migration is driven mostly by 
internal perturbers (Appendix, Sec. 19.11) . We show in 
the Appendix (Sec. 19.21 Fig. [T3")) that a sizable fraction 
of metal-poor, high-[0/Fe] stars can, nevertheless, mi- 
grate to the solar neighborhood as the effect of the bar. 
However, this happens over the entire evolution of the 
disc and the ratio of cold to hot number of objects is 
now ~ 0.9 (to be compared to 0.3 in our stander model), 
related to the overall low final velocity dispersions. In 
a properly mass-corrected observational data this ratio 
can be used as a discriminant for the two type of scenar- 
ios. A possible deficiency in the no-merger case is that 
velocity dispersion for the oldest stars is underestimated 
by a factor of ~ 2. Th is is consistent with the results by 
iMinchev et al.l (|2012b[ ) , in showing that internally driven 
migration does not contribute to disc heating. 

7. A UNIFYING MODEL FOR THE MILKY WAY THICK 

DISC 

We now focus our attention on the properties of the 
oldest stars (age> 10 Gyr) in our model. This subsample 
is marked by 

• a metallicity distribution which peaks at [Fe/H]~ 
—0.5 and has a metal-poor tail down to [Fe/H]~ 
— 1.3 (Fig. [J upper right panel). 

• [0/Fe]-values spanning the range 0.2 — 0.4, with a 
peak around 0.3 (Fig 0J lower right panel). 

• a lag in the rotational velocity by ~50 km/s com- 
pared to the young stars (Fig.[6j top row). 

• large velocity dispersions (Fig. [6]). 

• a large scale-height (Fig. bottom row) . 

All of these properties are strikingly reminiscent of 
what we call the "thick disc" of our Galaxy, despite the 
fact that we have used pure thin-disc chemistry. Within 
the frameworks of the model we present here, the MW 
thick disc has emerged from (i) stars born hot and heated 
by mergers at early times and (ii) radial migration (from 
mergers at early times and bar/spirals later on) trans- 
porting these old stars from the inner disc to the solar 
vicinity. 

Such a conclusion is in agreement with most (seemingly 
contradicting) models of thick-disc formation, which ex- 
pect contribution from only /mostly one of the follow- 
ing: (i) mergers, (ii) early formation in gas-rich, turbu- 
lent clumpy discs, or gas rich mergers, and (iii) radial 



migration driven by internal instabilities. A combina- 
tion of these mechanisms working together is required, 
where strong heating and migration occurs early on form 
external perturbations (our case) and/or turbulent gas 
clumps, followed by radial migration taking over the disc 
dynamics at later times. Yes, mergers are important, 
but we also need radial migration (unavoidable if a bar, 
spiral structure and/or mergers are present) to transport 
out old, hot stars, with thick-disc chemical characteris- 
tics. Yes, migration is important, but the old stars need 
to be "preheated" by being born hot and/or were heated 
by mergers at high redshift (also unavoidable from our 
current understanding of cosmology). 

The high stellar birth velocity dispersions at high red- 
shift we find in ou r simulation (~50 km/s) is consisten t 
with recent works ([Forbes et al.ll20ialBrook et al.l2012t ). 
An important dynamical consequence of this is that the 
disc becomes less susceptible to satellite perturbations 
(common at high redshift), making it easier to survive 
until today. 

Radial gradients for different height above the disc 
plane, as well as the scale-length of this "thick disc" pop- 
ulation will be discussed in detail in Paper II. 

8. CONCLUSIONS 

In this work we presented a new approach for study- 
ing the chemodynamical evolution of galactic discs, with 
special emphasis on the Milky Way (MW). Unlike other 
similar studies, where either the dynamics is too simplis- 
tic or the star formation history (SFH) and chemistry en- 
richment are unconstrained, our chemodynamical model 
is a fusion between a pure chemical evolution model and 
a high-resolution simulation in the cosmological context. 
As we argued in SeclU this new approach allows us to by- 
pass most known problems encountered in self-consistent 
simulations, where chemical enrichment still proves to be 
a challenge (see Sec. 11.21) . Moreover, this is the first time 
that a chemodynamical model has the extra constraint 
of defining a realistic solar vicinity also in terms of dy- 
namics (see Sec.[2|). 

The main results of our chemodynamical model can be 
summarized as follows: 

• The distribution of birth radii, ro, of stars ending 
up in a solar neighborhood-like location after 11.2 Gyr 
of evolution peaks close to ro = 6 kpc due to radial mi- 
gration (left panel of Fig. [3]). The strongest effect from 
radial migration is found for the oldest stars, related to 
the 1:5 mass-ratio merger in our simulation ~ 9 Gyr ago. 
Locally born stars of all ages can be found in the solar 
neighborhood. While a wide range of birth radii is seen 
for different age groups, the majority of the youngest 
stars are born at, or close to, the solar neighborhood 
bin. 

• While the low-end in our simulated metallicity dis- 
tribution is composed by stars with a wide range of birth 
radii, the tail at larger metallicities (0.25 <[Fe/H]< 0.6) 
results almost exclusively from stars with 3 < ro < 5 kpc. 
This is the region just inside the bar's CR, which is where 
the strongest outward radial migration occurs, as we dis- 
cussed in Sec. [2] and showed in the bottom row of Fig. [TJ 
The fraction of stars in this tail can, therefore, be related 
to the bar's dynamical properties, such as its strength, 
pattern speed and time evolution/formation. 

• Our results suggest that the most likely birth location 
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for the Sun is in the range 4.4 < ro < 7.7 kpc, with the 
highest probability ^5.6 kpc, followed by ^7 kpc (Fig.[3l 
right). This estimate comes form both dynamical and 
chemical constraints and is dependent on the migration 
efficiency in our simulation and the adopted chemical 
evolution model. 

• Examining the effect on the age-metallicity relation 
(AMR), we find that some flattening is observed, mostly 
for ages >5 Gyr (Fig. [4}. Although significant radial 
mixing is present, a slope in the AMR is preserved, with 
a scatter compatible with recent observational work. 

• We find no bimodality in the [Fe/H]-[0/Fe] stellar 
density distribution. However, when selecting particles 
according to kinematical criteria used in high-resolution 
samples to define thin and thick discs, we recover the ob- 
served discontinuity in the [0/Fc]-[Fe/H] plane (Fig. [2J) • 
This is i n agreemen t with t he recent observational re- 
sults by IBovv et al.1 (|2012bl ). where a smooth [Fe/H]- 
[O/Fe] distribution was obtained, after correcting for the 
spectroscopic sampling of stellar sub-populations in the 
SEGUE survey. 

• By separating our simulated local sample into narrow 
bins of age, [O/Fe], and [Fe/H], we found that the verti- 
cal scale-height of each component can be fitted well by 
a single exponential, with values growing with increasing 
age, increasing [O/Fe], and decreasing metallicity (Fig. 
bottom). The vertical velocity dispersions for each of 
these subpopulation was found to exhibit smooth varia- 
tions with height above the disc plane, strongly flattening 
for old, high- [O/Fe], and metal-poor stars. Our model re- 
produces qualitatively (and i n many cases quant itatively) 
the observational results bv lBovv et al.l ()2012al lcT). 

• We predict approximately fiat \z\ — [Fe/H] and 
\v z \ — [Fe/H] (perhaps anticipated) relations for stars of 
narrow age- or [0/Fe]-bins (Fig. [Sj) - This is an extension 
to the similar behavior seen in the |z| — a z ^ r plots above 
and provides additional predictions for Galactic surveys. 

• The local age-velocity relation (AVR) in the radial 
direction, ay, resulting from our model plateaus for stars 
of age 3-8 Gyr, while a weak increase is found for the 
vertical component, a z (Fig. [6|). For ages greater than 
8 Gyr, a continuous increase is observed for both com- 
ponents, related to the massive merger in our simulation 
~ 9 Gyr ago. These hot stars can be associated with a 
thick disc. In contrast, internally driven migration does 
not result in high enough velocity dispersions for the old 
stars (see Sec. EU Fig. H2|). 

• We found a strong flattening in the [Fe/H] radial 
profiles of the older populations, however, the younger 
ones are much less affected (see Fig. [5]). Remarkably, for 
stars younger than 2 Gyr the final gradient is very similar 
to the initial one out to ~ 12 kpc, justifying its use as a 
constraint for our chemical model. 

• We predict that the [O/Fe] radial profiles are es- 
sentially preserved for the chemical model we use. The 
[O /Fe] profiles for different age groups result straightfor- 
wardly from the adopted variation of the infall-law with 
radius (and hence the SFHs at different positions) and, 
thus, provide a way to constrain different chemical evolu- 
tion models. In the near future, these would be possible 
to measure by combining the good distances and age s 
expected from the CoRoT mission (BaghneFaD|2006), 
with abundance ratios obtained by spectroscopic follow- 
up surveys. For the young populations, this should be al- 



ready possible to be obtain from the observations of open 
clusters, e.g., with the ongoing Gaia-ESO or APOGEE 
surveys. 

• Finally, probably one of the most important out- 
comes of our chemodynamical model is that, although 
we used only a thin-disc chemical evolution model, the 
oldest stars that are now in the solar vicinity show sev- 
eral of the properties usually attributed to the Galactic 
thick disc. In other words, according to the results of 
the present work, the MW "thick disc" emerges natu- 
rally from stars migrating from the inner disc very early 
on due to strong merger activity in the first couple of 
Gyr of disc formation, followed by further radial migra- 
tion driven by the bar and spirals at later times (see 
Sec. E3). 

We showed that hot, old, metal-poor, high-[a/Fe] stars 
can be delivered to the solar neighborhood at high red- 
shift as the effect of a massive perturber (1:5 mass ra- 
tio in our simulation, see Fig. [TTj) . For a MW-size 
galaxy, a merger of this size around z ~ 1 is consis- 
tent with the numerical results by a number of groups 
(e.g.. IBenson et al.l 2004; Stochr 200 61: IKazantzidis et all 
| 2008t IDe Lucia and Helmil 120081: IStewart et all 120081: 
IVillalobos and Helmill2008t IPurcell et al.ll2009Q . 

Alternatively, a bar can also bring some fraction of 
these stars to the solar neighborhood, but over a large 
period of time, as we showed in Sec. l9.2l in the Appendix. 
However, a difference in the metallicity and velocity dis- 
tributions of these old stars would result. Most notably, 
avoiding a merger not only avoids the extreme migration 
at high redshift, but also does not result in high enough 
velocity dispersions today. For example, we find AVRs 
with maxima at oy ~ 43 km/s and a z ps 25 km/s in 
the lack of strong mergers (to be compared to ps 60 and 
ps 50 km/s, respectively). Additionally, the ratio of cold 
to hot orbits changes from ~0.3 to ~ 0.9 when the merger 
is avoided. This larger difference can be used as a con- 
straint in properly mass-corrected observational data to 
assess the possibility of an early-on merger in the MW. 

It is, therefore, tempting to conclude that stars with 
thick-disc chemistry and kinematics have migrated to the 
solar vicinity early on as the result of a sufficiently mas- 
sive merger. For an efficie nt migration, we sho uld expect 
a merger on a radial orbit (jQuillen et al.|[2009l ) penetrat- 
ing deep into the disc. Equivalcntly, it may be possi- 
ble that a series of sufficiently massive minor mergers 
(e.g., 1:10 mass ratio) acting over a short time period 
(~ 1 — 2 Gyr) may give a similar result. 

The bar in our simulations has a size at the final time, 
similar to the one seen in the MW, however, it is quite 
weak before the last ~5-6 Gyr of evolution and it does not 
show buckling, although the bulge appears boxy. It may 
be that our Galaxy hosted a stronger bar than the one 
we consider here, which would result in an even stronger 
migration from the inner disc. It would be interesting 
to explore disc models with different bar sizes, as well as 
cases where a strong buckling instability takes place. 

While here we concentrated on the solar vicinity, re- 
sults for the entire disc will be presented in Paper II 
of this series. Chemodynamical predictions for different 
disc radii may help further constraints thick-disc forma- 
tion scenarios and the MW assembly in general. 

The novel method presented here is potentially very 
useful and easy to implement. A suite of different chem- 
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ical evolution models assigned to different disc dynam- 
ics (from different simulations) should be investigated to 
make progress in the field of Galactic Archeology. 
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9. APPENDIX 
9.1. The effect of different model realizations 

As described in Sec. [2 we have downscaled our simulated disc by a factor of f=1.67, so as to place the bar's OLR 
just inside the solar circl e, at ~ 7.5 kpc (assuming the Sun is at 8 kpc), in agreement with its expected effect on the 
local disc dynamics (e.g., [Dehnen 2000; Minchev et al.H20"Tof ). 

We now would like to see how our results change if we rescale the disc by a different factor. We use f=1.4, which 
places the solar radius ~ 1.25 kpc inwards of the standard location used in the paper, which is equivalent to having a 
stronger (longer). We refer to this realization as Model A. 

As a second test, we start the implementation with chemistry 2.7 Gyr later, thus avoiding the early massive merger. 
We integrate the simulation for additional 2.7 Gyr so that we again have 11.2 Gyr of evolution. To keep the correct 
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location with respect to the bar's resonances, we downscale the disc radius by f=2.1 to account for the bar's slowing 
down. This is referred to as Model B. 

Note that in both cases above we only change the dynamics but not the input chemistry, i.e., we resample the SFH 
according to our chemical model, with the 8 kpc chemistry assigned to the rescaled 8 kpc radius in the simulation, etc. 

The results are shown in Fig. 1121 where in the first two rows we compare the birth-radius and metallicity distributions 
resulting from our new realizations with f=1.4 (Model A, first column) and the time shift (Model B, right column), 
with our standard rescaling (f=1.67, middle column). As expected, when the solar radius is shifted closer to the bar 
(Model A), a larger fraction of stars of all ages arrives from inside the bar's CR (dotted- red vertical line in top row); 
to see this more easily, we have overlaid the original total histogram (in pink) on top of the distributions. Inspecting 
the metallicity distributions, we see a shift in the peak of 0.1 dex to negative values when the Sun is closer to the bar. 
However, the changes are not drastic and the rest of the results in the paper are not affected much by this change. 
Whether this is true for the entire disc will be investigated in Paper II. 

As it became evident in the discussion of the age- metallicity relation f Sec. 15.3]) . in our simulation we have a deficiency 
of stars with ages>~ 9.5 Gyr at r > 10 kpc. However, this is a very small fraction of particles, as can be inferred 
from the bottom right panel of Fig. [2] (the SFRs at those times and radii attain a maximum at ~ 2 M© pc -2 Gyr -1 
for r — 10 kpc, sharply decreasing for older stars and larger radii). By shifting the solar radius inward, this artifact is 
now lost. How small this fraction is can be seen by comparing the extent of the red ro-histograms (oldest stars) in the 
first two rows outside the solar bin. 

The right column of Fig. [12] shows the new realization with the time shift (Model B), where we have avoided the 
early-on massive merger, but keep the bar resonances the same distances from the Sun. A smaller number of stars 
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Fig. 13. — Comparison between the effect of an early massive merger (first two columns) and internally-driven radial migration (third 
and fourth columns) on the dynamics of the oldest stars. See text. 

arrives from inside the bar's CR now, compared to our standard case (see overlaid pink line). This comes mostly from 
the lack of the strong peak in the ro-histogram of the oldest stars (red). We show below (Fig.[T"3j) that the enhancement 
in the peak at ro ~ 3 kpc is due to the merger, which we avoid here. 

In the bottom right panel of Fig. Q2] we compare the vertical velocity dispersion, a z (r) (solid solid), and rotational 
velocity, V^(r) (blue dashed), radial profiles resulting from the new realizations with our standard model (pink curves). 
A slightly larger velocity dispersion is seen for the inner disc radius (left), while the difference is drastic when the 
massive merger is avoided. Most importantly, the oldest stars have o~ z about a factor of two smaller than when the 
merger is present. The maximum radial velocity dispersion values (not shown) also drop to ~ 43 km/s, which is 
reflected in the higher V^-values. These smaller velocity dispersions and rotational velocity difference between young 
and old stars is unlike the observed values in the solar neighborhood. 

9.2. The effect on the dynamics of the oldest stars in the absence of a massive early-on merger 

In Sec. 16.31 we discussed the origin of the old, metal-poor, high-[0/Fe] stars in our simulated solar neighborhood. 
We showed in Fig. [TT] that these stars are born in the range 1 < r < 10 kpc, where a large fraction is already in place 
at t = 3 Gyr, i.e., shortly after the last massive (1:5 mass ratio) merger ~ 9 Gyr ago. 

We now examine the fate of these stars in the absence of a massive early-on merger (Model B in Sec. I9.1j) . The first 
two columns of Fig. 1131 show the birth radius distributions, ro, for the hot (first row) and cold (second row) orbits at 
the final state (top) and at t — 3 Gyr (bottom) for the model realization used in the paper. These are the same as 
the bottom two rows of Fig. 1111 In the third and fourth columns we now plot the equivalent subsamples, but for the 
model realization lacking any massive mergers (Model B). We first compare the final states (top row). The almost 
square shape of the distribution for the hot orbits on the left is now replaced by a triangular shape with a maximum 
just inside the bar's OLR, i.e., less old stars from the inner disc end up near the Sun when the merger is avoided. 

As we showed in Fig.[T2J bottom right, for Model B the oldest stellar population has velocity dispersions lower than 
the standard case. Therefore, separating stars by a cut in eccentricity at e = 0.2 results in a "hot" sample with a 
radial velocity dispersion oy sa 55 km/s, to be compared to « 65 km/s (on the other hand, both cold samples have 
oy r; 28 km/s). The ratio of cold to hot orbits changes from ~0.3 to ~ 0.9. This large difference can be used as a 
constraint in observations to assess the possibility of an early-on merger in the MW. 

Concentrating on the bottom row of Fig. [13] we now see that although it starts with a fully formed bar, Model B 
is unable to mix the disc by t = 3 Gyr as much as in the case of the early-on merger. Only a small fraction of stars 
arrives from inside the bar's CR for both the hot and cold orbits. 

In conclusion, internally-driven radial migration (compared to the case of an early-on merger) results in (i) a much 
smaller fraction of old stars originating from r < 4 kpc and (ii) these stars are transferred continuously throughout 
the disc evolution. 

9.3. Migration contribution to the disc velocity dispersions during and after the massive merger 

We here shed some light on the effect of migration on the changes in velocity dispersion radial profiles during and 
after the early-on massive merger in our standard model. 

We separate migrating from non-migrating stars (in a given period of time) by applying the technique described 
bv lMinchev et al.l (|2012bD . Fig. [T4l shows the migrators' contribution to the disc heating during the merger (1 < t < 
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1<t<2.5 Gyr 




Fig. 14. — Migrators' contribution to the disc heating during the time period the massive merger takes place (first and second columns) 
and during a period immediately after that (third and fourth columns). First and Third columns: Radial profiles of the vertical (top) 
and radial (bottom) velocity dispersions, cr z (r) and <x r (r), for non-migrators, outward going migrators, and inward migrators (as indicated) 
for each time period. Second and Fourth columns: The effect of non-migrators, all migrators (i.e., inward and outward going), and 
the total population. Note that during the merger (left two columns) velocity dispersions are higher for stars coming from the outer disc, 
followed by the non-migrators, while samples originating in the inner disc are the coolest. This is due to the stronger merger effect on the 
outer disc. In contrast, at later times (right two columns), the bar and spirals transfer the already preheated stellar populations, thus, 
contributing to disc heating. 

2.5 Gyr, first and second columns) and during the time period immediately after that (2.5 < t < 4 Gyr, third and 
fourth columns). The radial profiles of the vertical (top) and radial (bottom) velocity dispersions, cr z (r) and <7 r (r), for 
non-migrators (dashed black curve), outward going migrators (solid brown curve), and inward migrators (dash-dotted 
green curve) for each time period are shown in the first and third columns. The second and fourth columns show the 
effect of non-migrators (dashed black curve), all migrators (i.e., inward and outward going, solid green curve), and the 
total population (dashed-dotted blue curve). 

During the merger (left two columns) velocity dispersions are higher for stars coming from the outer disc, followed 
by the non-migrators, while samples originating in the inner disc are the coolest. This behavior (especially for a z ) 
is opposite to the effect expected from internally-driven migration (Mi nchev et al.ll2012bf) and is due to the stronger 
perturbation experienced by the outer disc, given its exponential surface density decrease. The second column shows 
that, for r > 3 kpc, the overall effect of migrators is to decrease the total velocity dispersion, because most migrators 
come from the inner disk. 

In contrast, at later times (right two columns), migration induced by the bar and spiral structure transfers out the 
already preheated stellar populations, thus, contributing to disc heating. As the disc grows further from gas accretion, 
this contribution becomes much smaller. The weak effe ct of migration o n disc heating at the later stages of this same 
galaxy evolution was presented in Fig. 8 (model C3) by lMinchev et al.l (|2012ri ). 



